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, I Abstract. A mass ejection model in a time-dependent random environment with 

1*^ both temporal and spatial correlations is introduced. When the environment has a 

r — . finite correlation length, individual particle trajectories are found to diffuse at large 

04 times with a displacement distribution that approaches a Gaussian. The collective 

^_^ dynamics of diffusing particles reaches a statistically stationary state, which is 

^^ characterized in terms of a fiuctuating mass density field. The probability distribution 

r \ of density is studied numerically for both smooth and non-smooth scale-invariant 

random environments. A competition between trapping in the regions where the 
ejection rate of the environment vanishes and mixing due to its temporal dependence 
leads to large fiuctuations of mass. These mechanisms are found to result in the 
presence of intermediate power-law tails in the probability distribution of the mass 
density. For spatially differentiable environments, the exponent of the right tail is 
shown to be universal and equal to —3/2. However, at small values, it is found to 
depend on the environment. Finally, spatial scaling properties of the mass distribution 
are investigated. The distribution of the coarse-grained density is shown to posses 
l/^ some rescaling properties that depend on the scale, the amplitude of the ejection rate, 

fK^ and the Holder exponent of the environment. 
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1. Introduction 

Many different situations found in nature take place in an environment where spatial 
fluctuations occur at so small length scales that they contribute only through an 
averaged effect on macroscopic processes. The environment is then usually modeled 
by introducing some disorder, like in the Sherrington and Kirkpatrick spin-glass model 
for magnetization [I]. In addition the time variations of the environment are often 
much slower than those of the mechanism of interest (such as diffusion, transport or 
wave propagation). This leads on to consider a quenched disorder, as in the case of 
directed polymers [21 E] and of wave propagation in random media [1]. As a consequence, 
the discrete version of diffusive processes in quenched random media, known as random 
walks in random environment (RWRE) [5|, have been largely studied by mathematicians 
and physicists. One usually considers a lattice with fixed random transition probabilities 
between sites and study the behavior of a random walk on it. One of the main 
questions of interest has been to determine under which assumptions such RWRE 
are transient or recurrent (that is whether they escape to infinity or indefinitely come 
back to their starting position). For a time-independent random environment with no 
space correlations, it has been shown that randomizing the environment slowdowns 
the diffusive properties of the random walk [6]. Furthermore, under some precise 
assumptions, Ya. Sinai [7] proved that one-dimensional symmetric nearest-neighbor 
random walks in an uncorrelated environment are sub-diffusive and scale as Xt ~ log^ t 
when t — i- oo. B. Derrida and Y. Pomeau [S] demonstrated that when the assumption 
of symmetry is relaxed, the random walk escapes as Xt ~ t", with < a < 1. 
J. Bricmont and A. Kupiainen [0] showed that the upper critical dimension of random 
walks in a time-independent random environment is 2, so that ioi d > 2 RWRE are 
diffusive. More recently, a large effort has been devoted to prove central-limit theorems 
or large-deviations principles for quenched and annealed disorders [ini [HI [12] • Such 
work generally assume that the environment has no (or very short) spatial correlations 
and that it is time- in dependent. Only little is known for the case of time-dependent 
environments. In such settings and under some rather general assumptions, it was shown 
that random walks are always diffusive [13] and that central limit theorems generally 
hold [H]. 

We are here interested in situations where the timescales and length scales of 
diffusion are comparable or longer than those at which the environment fluctuates. 
A clear instance where this occurs is turbulent transport whose modeling is a key issue 
in industrial and environmental sciences. The classical models used for instance in 
engineering and meteorology are based on an "eddy-diffusivity" approach (see, e.g., [15]). 
which in general requires a clear scale separation between the turbulence and large-scale 
variations of the mean velocity. Homogenization techniques can then be used to show 
that the averaged concentration field follows an effective advection-diffusion equation 
with an advecting velocity and a diffusion coefficient tensor that depend on the slow 
variables [16]. The large-scale mixing originating from small-scale fluctuations acts only 
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through the diffusive term, which, by the maximum principle, cannot be responsible for 
the creation of large concentrations. The concentrations observed in compressible flows 
can thus only come from the effective advection term. As we will see in this paper, 
the situation can be rather different when interested in other types of compressible 
transport. 

Let us consider for instance species whose dynamics is dissipative in the full position- 
velocity phase space. We assume that the trajectory of a particle obeys the Newton 
equation 

Xt = Vt, Vt = -fiVt + ^fiXt,t,t/6). (1) 

The forces acting on the particle are a viscous drag and an external force / that depends 
on both space and two different time scales. We have dropped here the vectorial 
notations but all the following considerations can be easily generalized to any dimension. 
Such an equation describes for instance the dynamics of a heavy inertial particle in a 
velocity field that varies over two time scales. Also, one could consider a large particle 
embedded in a time and space-dependent thermal bath. The fast time scale e can be 
interpreted as the typical time of momentum exchange between the particle and its 
environment. As in the Einstein original work about Brownian motion, the time-scale 
e is much smaller than the viscous damping time l//i. In the limit e/i — )■ 0, the force / 
can be approximated by a Gaussian noise with correlation 

C(t, t') = (l/e) (/(X„ t, t/e)fiX,,t\ t'/e)), ^ 2Tfi'a\X,, t)6it - t'), (2) 

where the average {■)s is with respect to the fast time variable. The spatio-temporal 
variations of cr^ can be interpreted as a non-homogenous temperature field in the thermal 
bath. We next consider the limit when the response time 1/fi is much shorter than the 
slow time scale. This introduces a new fast time scale 0(l//i) over which the particle 
velocity fluctuates. We see that 

t t' 

{VtVt,),= I f C{s,s')e^'^'+''-'-''Msds' 



— oo — oo 
tAt' 



(3) 



where the last equality is obtained integrating by parts and neglecting exponentially 
small terms. Then, taking the limit /i — )■ oo in pi) we obtain (V* Vf) = 2Ta^(Xt, t)5{t — 
t'). The particle position thus satisfies the stochastic equation 



dXt = V2ra{Xt,t)dWt, (4) 

where Wt is the Wiener process and where the product has to be interpreted in the Ito 
sense. The averaged particle concentration p is then a solution to the associated forward 
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Kolinogorov equation, namely 

|-r|,[.^,M)p]. (5) 

The solutions to such an equation have a very different behavior than those to the 
standard diffusion equation where on the right-hand side a'^ would appear in between 
the two spatial derivatives. It is for instance clear that the mass is going to accumulate 
at the zeros of cr. Indeed, suppose a{x,t) ~ Cx in the vicinity of x = 0. At leading 
order, the flux reads J = —T{d/dx){a^p) ~ —2TC'^xp, which is positive for a; < and 
negative for x > 0, leading to a permanent mass flux toward x = 0. As we will see 
later, the zeros, their densities and their lifetimes play a crucial role in the statistical 
properties of the density fleld. 

The paper is organized as follows: in section |2] we introduce a time and space 
continuous model of diffusion in a random environment based on an ejection model in a 
discrete D-dimensional lattice. The model is reinterpreted in terms of particle dynamics 
obeying an Ito differential equation with a multiplicative noise and no drift. The random 
environment is then explicitly deflned and its general spatial and temporal properties 
are discussed. Section |3] is devoted to the study of individual trajectories of the Ito 
diffusion in the case of smooth and non-smooth environments. When the environment 
has a flnite correlation length, particles are found to have standard diffusion properties 
at large times with a displacement distribution that approaches a Gaussian. In section 
|4]we deflne the density of mass and study the fluctuations of the density in smooth and 
non-smooth random environments. The locations where the environment ejection rate 
vanishes are shown to play a crucial role on the large fluctuations of density. Spatial 
scaling properties of the mass distribution are then investigated in [5] We show that 
the distribution of the coarse-grained density possesses some rescaling properties that 
depend on the scale, the amplitude of the ejection rate, and the Holder exponent of the 
environment. Finally, Section |6] presents some concluding remarks. 

2. Description of the model 

Let us consider the D-dimensional periodical lattice of fundamental size Ax and period 
L = NAx. This lattice deflnes a tiling on which we consider the following discrete-time 
dynamics. The cell indexed by i contains at time nAt a mass Pi{n). Then, between 
times nAt and {n + l)At, a fraction < 7i(n) < 1/{2D) of this mass is ejected from the 
cell i to one of its 2D neighbors. For a given set of ejection rates {7i('^)}ie{i,...,Ar}-o the 
variables {pi(^)}i6{i,...,Ar}i3 deflne a Markov chain whose master equation is given by 

Pi{n + 1) = [1 - 2D^,{n)] p,{n) + ^ 7j(^)Pj(^), (6) 

jeM 

where A/l are the neighbor cells of i. It is clear that, because of periodic boundary 
conditions, the total mass M = J2i Pi is conserved. Equation ^ describes a mass- 
ejection process in a time-dependent non-homogeneous environment determined by the 
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ejection rates 'yi{n). A similar model has been studied in [17] in the case when the 
7i(n) take two values, 7 and and are un correlated in both space and time. This 
model was motivated by the study of inertial particles clustering in turbulent flows 
and was expected to mimic heavy particle ejection from coherent vortical structures by 
centrifugal forces. 

Turning back to the ejection model ([6]), we now consider the continuous limit 
Ax — )■ and At — t- 0. To take this limit, we suppose that the ratio Ax^/At is kept 
fixed. We denote by ro"^(x, t) the continuous limit of 'y:x_/Ax(t/At), where the coefficient 
r is the typical amplitude of the ejection rate and a is an order-one time and space 
continuous function. In this limit, equation ^ becomes 

^^ = rV^[a^(x,t)p(x,t)]. (7) 

As mentioned in the introduction, in contrast with a standard diffusion equation coming 
from Pick's first law, the differential operator V^[cr^(x, t) ■] in the right-hand side of 
equation (ul) is not positive-definite. There is no maximum principle and the solution is 
expected to behave very differently from those to standard diffusion equations. 

Equation (JTl) is the Fokker-Planck (forward Kolmogorov) equation associated to the 
Ito stochastic differential equation 



dXt = V2ra{Xt,t)dWt, (8) 

where Wt is the D-dimensional Wiener process. For t > s, we define the (forward) 
transition probability density 

p(x,t|y,s) = (5(Xt-x)|X, = y), (9) 

where (•) designates averages with respect to the realizations of W^. Then the solutions 
to ([7]) trivially satisfy for allt > s 

p(x, t) = p(x, t\y,s) p(y, s)dy. (10) 

The model described by the diffusion equation ([T]) can thus be reinterpreted in terms of 
a random walk in a random environment (RWRE). Indeed the discrete version of the 
process given by (|8| corresponds to a random walk on the D-dimensional lattice, with 
a probability to hop from the site i to one of its neighbors equal to 2D'ji{n). 

Let us now specify the statistical properties of the random environment. We assume 
that L = 271 and that cr(x, t) is 27r-periodic in space. The environment is then entirely 
determined in terms of its Fourier series that we write as 

a(x,t) = 5^akXk(t/rk)e^'^^ (11) 

|k| 

where the a^'s are positive real amplitudes and the Tk's are scale-dependent 
characteristic times. The Fourier modes Xk satisfy X-k(^) = Xk(^) ^^^ Xo = 0. They are 
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independent Gaussian processes with unit variance and unit correlation time. We choose 
them as complex Ornstein-Uhlenbeck processes that solve the stochastic differential 
equation [IB] 

dXk(s) = -Xkds + v^dfik(s), (12) 

where the -Bk's are independent one-dimensional complex Wiener processes. 

We now prescribe scale- invariance properties for the random environment. We want 
|(j(x + r) — o-(x)| ~ |r|'', where h > controls the regularity of the ejection rate. When 
/i < 1, it is the spatial Holder exponent of the ejection rate. This amounts to assuming 
that the Fourier mode amplitudes behave as a^ = ^fc~^/^~'*. Temporal scale invariance 
is set by assuming r^ = k~^ with /3 > 0. The exponent /3 relates time dependence to 
spatial scale invariance; the case (3 = corresponds to a having a unique correlation 
time and /3 — )■ oo to white-noise. As we shall see below, when h < 1, dimensional 
analysis motivates the choice /3 = 2 — 2h. 

The numerical results presented in this work are obtained by two different manners. 
The density evolution is obtained by solving equation ([t]) with a second-order finite- 
difference discretization of the Laplacian and a semi-implicit Euler temporal scheme. 
This numerical scheme ensures the conservation of the total mass with a high accuracy. 
In this work numerical resolutions vary from N = 512 to A^ = 8192 collocation points 
and the time-step is chosen small enough to well resolve the fastest time scale of the 



problem. For particle trajectories, the stochastic equations (\8^ and (12) are solved by a 



standard stochastic Euler scheme and ensemble averages are obtained by Monte-Carlo 
methods. In both cases, we expect the error made on the solution to act as a numerical 
diffusion with a constant proportional to the time step. 



Density q 



Ejection rate o 





{a)D 



(h) D = 2 



Figure 1. (a) Left: space-time plot of p{x,t) (liigli densities in yellow and low in 
black) for /i = 1 in dimension D = I a Right: corresponding evolution of the ejection 
rate a^{x,t) during the same realization; low values are in dark and high values in 
yellow, (b) Snapshot of the position of particles obeying ((sl) (in green) for _D = 2, 
together with a'^{x,t). Low values of a'^{x,t) are in dark and high values in yellow; 
contour lines of cr = are shown in white. 
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Figure |l(a)| show the temporal evolution of the density and of the ejection rate 



cr^(x,t), respectively, for a typical realization in dimension D = 1. Figure 1(b) shows a 
snapshot of the random environment a and of the position of a set of particles obeying 
equation ([s]) in D = 2. It is apparent on figure [I] that the high density zones are located 



in the vicinity of the zeros of o"(x, t). From the right of figure 1(a) we observe that the 
zeros of a follow random paths. Their role on the statistical properties of the density 
of mass and on the diffusion properties of particles will be further discussed trough this 
work. 

The distribution and time-evolution of these zeros strongly depend on the parameter 



h characterizing the environment. Let us first consider the case of large h. On figure 1(a) 
it is apparent that the zeros typically appear by pairs and then diffuse and separate until 
they merge together or with another zero. Without loss of generality, we can assume 
that the ejection rate has only two modes and takes the simple form 

a{x,t) = (Jr(t) cosx — (Tj(t) sinx, (13) 

with o"r + icTj = 2%!. Using the Ito formula it is possible to show that a can be 
equivalently rewritten as a{x,t) = At sin (x — $j), where the stochastic processes At 
and $i satisfy 

dAt = -Ut- ^) dt + V2dBf , d$, = ^di?f , (14) 

where Bf^ and Bf are uncorrelated Wiener processes. It is clear that the amplitude 
At fluctuates around A = 1 and that its correlation time is order one. It follows from 



equation (14) that, for typical values of At, the phase $t, and thus the position of the 
zeros of a, diffuse on a order-one timescale. 

When h < 1, the random environment presents scaling properties and a very 
different behavior is expected. We have by construction that the second-order structure 
function of a behaves as 



Saj = \a{x,t) - a{x + £,t)\^ ^ {i/Ly\ (15) 

where the over-line stands for the ensemble average with respect to the fluctuations of the 
environment (i.e. with respect to the Ornstein-Uhlenbeck processes Xk)- In this work, we 
make the choice of relating space and time correlations by a dimensional argument. The 
characteristic time r^ introduced in equation ( [ll| ), behaves as a power law, so that the 
correlation time of a at scale £ is tc{£) = {£/ L)^ . According to equation ( 15 ) the diffusion 



timescale associated to the spatial scale £ behaves as td{£) = £'^ /{TSaf) ~ T^^ L"^^ £'^^'^^ . 
The ratio Ku(^) = tc{£^) / tj^{£) ~ p/^-(2/i+/3)£/3+2/i-2 (;jggj^gg ^ dimensionless space- 
dependent parameter that is usually called the Kubo number. When Ku(£) ^ 1, 
the environment is as if frozen. When Ku(£) ^ 1, it fluctuates in an almost time- 
uncorrelated manner. When /3 7^ 2 — 2/i, the Kubo number depends on the scale and 
this breaks any possible scale invariance of the mass distribution. When /3 = 2 — 2/i, 
we have Ku(£) = Ku = T/L^, so that scale invariance is possible. Here, we focus 
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on this latter case and will describe the mass concentration properties in terms of the 
dimensionless parameter Ku oc F. Note that the choice of having a scale-independent 
Kubo number is common in the framework of turbulence and is expected to be relevant 
to the problem of diffusion of inertial particles in turbulent flows. 

Finally, let us do a couple of remarks on the dependence upon the other parameter 
of the environment, namely the Holder exponent h. Using Parseval's theorem, it is 
possible to show that ||cr||^2 = C(l + 2/i) and that ||92:Cr||^2 = C(^l + 2/i), where ({s) 
is the Riemann zeta function. We thus have that ||cr||/^2 < oo when h > and that 
||<9a;C'"||i2 < CO when h > 1. Note that, at each time t fixed, a{x,t) is a Gaussian 



variable with the variance of increments given by equation (15). Therefore (j{-,t) is 
a fractional Brownian motion of exponent h fl9\. For h = 1/2 it corresponds to the 
standard Brownian motion. For h > 1/2 the increments of a are not independent and 
have a positive covariance; the zeros of a are then finite and isolated. For h < 1/2 
the covariance of increments is negative and therefore we expect the number of zeros to 
become infinite and to accumulate. This different behavior of the random environment 
will affect the general properties of the diffusion, as we will see in the next sections. 

As the Ornstein-Uhlenbeck processes Xk are stationary, the fluctuations of the 
ejection rate a"^ is a stationary and homogeneous random field. Hence we expect that, 
at sufficiently large times, the density distribution reaches a statistically stationary state 
and numerical simulations indicate that this is indeed the case. The work reported in 
this paper mainly concerns with the statistical properties of the density field in this 
large-time asymptotics. Most of the results presented in this work are in dimension 
D = l. 

3. Diffusive properties 

We now turn to studying the diffusive properties of the solutions X{t) to ^. One has 

X{t) = X(0) + v^ f (t{X{s), s) dW„ (16) 

Jo 

so that 

([X(t)-X(0)]^>=2r f\a'iXis),s))ds, (17) 

Jo 

where the angular brackets denote average over trajectories. At times much larger than 

the correlation time of a^ along particle trajectories, the integral becomes a sum of 

independent random variables, which obeys the law of large numbers. We thus have 

{[X{t)-Xm')-2D{T,h)t 

with D{T, h) = T{a^{X{t),t)) = \ f a^{x, t)p{x, t) dx, (18) 

-^ Jo 

where the over-line denotes here the average with respect to the environment. The 

effective diffusion constant -D(r, h) involves an average of the environment along particle 

trajectories which is equivalent to an average weighted by the mass density p. As the 
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trajectories are expected to spend a long time in the vicinity of the zeros of a, we 
expect this Lagrangian average to be less than the full Eulerian spatial average. The 
displacement fluctuations of order less than \/t are expected to be given by a central 
limit theorem. However, larger fluctuations should obey a large deviations principle 
with a rate function that might not be purely quadratic as it is non trivially related to 
the Lagrangian properties of the environment. To study these fluctuations, we perform 
Monte-Carlo simulations of equation ^. 

We first consider the case where the ejection rate a"^ is a smooth function of space 
(i.e. h > 1). With this configuration there are a finite (and small) number of isolated 
zeros separated by a distance order L/2 = n. A number of simulations of equation (Is]) 
have been performed for different values of F using 32 modes and h = 2. Figure |2] shows 
the mean-square displacement of particles averaged over 1000 realizations of the diffusion 



and 20 000 realizations of the environment. As expected from equation (18), a linear 
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Figure 2. Mean-square displacement {iX{t) — X{0))^) for different values of F. 



Inset: Dependence of the effective diffusion constant (18 1 on F 



growth is observed at large times. However the diffusion constant -D(F) non trivially 
depends on F, as seen in the inset of figure [2] The non-monotonic behavior of D{T)/T 
can be explained by two competing phenomena. In the limit F — )■ (or equivalently 
Ku — )■ 0) the environment changes fast enough compared to the characteristic time 
of diffusion given by 1/F. In this limit, we can consider that the particle trajectories 
sample homogeneously the environment, so that the Lagrangian average can be replaced 
by the spatial Eulerian average. This leads to 



^(F, h)-- a{x, tf dx = FC(1 + 2h) when 7 -> 0. 



(19) 



For h = 2, we obtain i5(F,2)/F ^ 1.03, which is in agreement with the value observed 
on the inset of figure |2} As F increases, trajectories spend a longer time at the zeros of a 
decreasing the value of D(F)/F. In the limit of F — )■ 00 (i.e. Ku — )■ 00), the environment 
can be considered completely frozen and we expect the mass density solving equation 
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Figure 3. Rescaled PDFs of the displacement X{t) — X{Q) for h — 2, for three 
different values of T and for various times, as labeled. The initial condition was chosen 
in the statistically stationary regime. 



([T]) to be approximatively given by p{x,t) ~ l/a^(x,t). From equation (18) we get 
^(r, h) = {T/L)J^dx ^ iV/L) / dx = r. Hence D{T)/T = 1 for T ^ oo as seen on 
the figure. In this hmit the mass is completely concentrated at the zeros of the ejection 
rate, and the diffusion is carried out by the motions of the laters that follow a dynamics 



similar to equation (14). 



The normalized probability distribution function (PDF) of the displacements are 
shown on figure |3] for different times and different values of F. At early times and 
for low values of F, the PDF's present exponential tails, that are also observed at 
intermediate times for larger values of F. In all cases, the PDF's approach a Gaussian 
distribution at very large times. Note that for sufficiently large values of F, some 
oscillations can be observed on the PDF tails at intermediate times. As we will now 
see, this is due to trapping by the zeros of cr(x, t). To emphasize this point, we compare 
the PDF of the displacement for different values of F chosen at the time t* such that 
(((X(t*) — X(0))^)) = {L/2Y = TT^, i.e. the time when the typical distance traveled by 
the particles reaches a length of the order of the separation between two zeros (of the 
order of L/2 for smooth environments). The PDF of {X{t*) — X{0))/{L/2) is displayed 



on figure 4(a) for various values of F. It is apparent that the bumps appear at multiples 



of L/2 for large F's. With a high probability, a particle is initially located close to a 
zero of a. When it travels away, there is a large chance that it stops in the neighboring 
zero, leading to a quasi multimodal distribution. 

To better quantify the large-time behavior of the displacement, we compute 
its higher-order moments and compare them to those corresponding to a Gaussian 
distribution. We thus define for even p 



Gri 



((x(t) - x{o)y) 2P/2 



((x(t)-x(o))2)p/2 y^ 



p + i 



(20) 
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Figure 4. (a) PDF of {X{t) - X{0))/{L/2) for different values of F taken at t* 
such that {{iX{t*) - X{0)f)) = iL/2f = tt^. (b) Temporal evolution of normalized 
moments Cp. defined by (|20]) for F = 0.1 (solid lines), F = 1 (dashed lines), and F = 10 



(dotted lines), (c) A{T) showing the dependence of Cp on F, see equation (21 1. 



where r[z] denotes here the Gamma function. Cp = for all p > 2 if X{t) — X(0) is a 
Gaussian variable. C2 = by construction and C4 is the Kurtosis. Cp can be interpreted 
as a p-th order deviation from a Gaussian distribution. Data suggest that all the Cp's 
display a 1/t behavior at large times. Its dependence on the amplitude F and on the 
order p is well represented by the functional form 



C ~ 



A{r)xp/^ 



(21) 



The temporal evolution of Cp/A{r) for different values of F and p, and the function 
A{T) are displayed on figures [4(b) and |4(c) respectively. As seen from 4(b), the data 
for different F and p fixed collapse at large times. The limiting lines associated to 
different orders are equally space, confirming the dependence in A^/^ of the constant. 
The value of A was found to be A ~ 15. 

Different simulations have been performed for various values of the exponent h. 
A behavior similar to the smooth case was found: at large times particles diffuse and 
their statistics become Gaussian. The rescaled PDF's of the displacement are shown 
on figure [5(a)] and 5(b) [for F = 1 and for h = 0.2 and h = 0.6, respectively. At large 



times {{X{t) — X(0))^) also present a linear grow consistent with equation (18) (data 



not shown). The measured diffusion constant D{T,h) is displayed on figure [5(cy for 
F = 1 and different values of h. The apparent divergence at small values of h is due to 



the divergence of the L^ norm of a{x,t) for /i — )■ 0. The inset of figure 4(c) shows the 
dependence of the diffusion constant normalized by Dq = ||cr||?2 = C(l + 2/i). It is seen 



that the importance of zeros on the diffusion constant increases when h decreases. 



4. Density fluctuations 

So far, we have studied the diffusion properties of particles described by the stochastic 
equation (JSJ). We now turn to studying the statistical properties of the mass density p 
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Figure 5. Rescaled PDF of X{t) - X{0) for F = 1 and with (a) h = 0.2 and (b) 
h = 0.6 for different times. Dashed hnes correspond to the Gaussian distribution, 
(c) Diffusion constant for F = 1 as a function of h. The inset show the same data 
normahzed by Dq — ||cr||^2- 

that evolves according to the diffusion-hke equation (JTl). As we have previously observed, 
the spatial fluctuations of p are correlated with the distribution of zeros of the random 
environment 0"(x, t). More precisely, we observe from equation (JTl) that the flux of mass 
at X is directly given by 



J(x, t) = r a(x, t) [2Va(x, t)p(x, t) + or(x, t) Vp(x, t)] 



(22) 



The flux hence vanishes at the zeros of cr(x, t), so that they act like sinks and concentrate 

all of mass if they were not moving. However, due to their diffusion, the zeros are not 

able to indefinitely concentrate mass and the density saturates in their neighborhood. In 

parallel, as the total mass is conserved, this concentration process leads to the creation 

of voids in the regions where a is order one. 

In the following we consider the case D = 1. Under some assumptions of 

homogeneity and ergodicity of the dynamics, the PDF of the mass density can be written 

as 

d 



P{p) 



dp 



(measured (x,t) \p{x,t) < p})) 



T 271 



lim — — 

T^oo 2'kT 




5 {p{x,t) — p) dxdt 



lim Prip) 



(23) 



(24) 







This will be the main quantity studied in the next sections. As we will see it strongly 



depends on the parameters F and h. The definition (24) in terms of a space-time 



average will be particularly useful to attack first the stationary case and then to develop 
phenomenological arguments on its behavior in the non-stationary case. 



4-1. Smooth random environments 

Let us first consider the limit Ku — )■ oo, i.e. when the ejection F is infinitely larger 
than the inverse of the diffusion time of the zeros of a. Note that, in the limit of a 
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time- independent environment, no stationary state can be achieved. Mass will then 
concentrate in the zeros and at t = oo the density will become atomic, supported on 
these points. For sufficiently large but finite times, as almost all the mass is concentrated 
around the zeros of a, the ejection rate can be Taylor-expanded in the vicinity of these 
points. Without loss of generahty, let us assume that a zero appears at the origin x = 
at t = and does not move. In this limit, equation ([T]) reduces to 



dtp = TC' 






[x^p), 



(25) 



where C = da/dx\x=o- Rescaling time allows us to set T C^ = 1. Without loss of 
generality we moreover consider an infinite domain and an initial condition with compact 
support: p{x,0) = 1 for |x| < 1 and elsewhere. Equation (25) can be analytically 



integrated by using the method of characteristics. For x > the change of variable 
u{y,t) = e^^p{x = e^~'^*,t) leads to the homogeneous heat equation that can be easily 
solved by using the corresponding Green function. The solution for a; < is obtained 
by symmetry. One then obtains 



p(x,t) 



where eTic{z) 



J2t 



'A-Kt 



exp 



{z — log |x| — 3t)^ 
4t 



dz 



3^*erfc 



In |x| -|- 3t 
2^^t 



(26) 



2 joo g ** ds is the complementary error function. The time-averaged 



V^ 



distribution Pt{p) defined in equation (24) can be rewritten as 



r 27r 



Pt{p) 



27iT 








S{x — Xp{t)) 

\ dp(x(t),t) I 
I dx I 



dxdt. 



(27) 



where Xpit) is such that p{xp(t),t) = p. It exists only for t > \Tip/2 . The contour 



lines of p{x,t) are displayed on figure 6(a) The high-density zones are concentrated in 
a narrow region of space. This will justify the saddle-point approximation used below. 
With the change of variables t i— )■ A = 2t/lnp, introducing p{X] p) = — In Xp{t)/ In p, 



and after some algebra, equation (27) becomes 



2T 
In p 



PT(p) = U'^Jv-Xe-'-'^"-^-^^MX, 



(28) 



with 



F(A;p) = A + /i(A;p) 
3 



1 
2A V 2 



/i(A;p) 



-A 



2A 



erfc 



-A-p(A;p) 



(29) 



^A-l 



;A-v/2A(A-1), 



Inp yp'" 

where in the last equation the asymptotic of eTic~^{z) ~ A/log (l/z) for z ^ has been 
used (recall that A > 1). Using a saddle-point approximation for Inp :^ 1 of the integral 



Mass fluctuations and diffusion in time- dependent random environments 



14 




10* 






















^ 


lo' 


" 


^ 


■< 




p ' 












" 


10° 










^^ 


-- 


-. 








- 


QlO"' 
Q. 


















X 


^, p"' 


- 


10* 




















Xn, 


~ 


10° 




















"^ 


^ 


10"' 






















- 



-0.15 -0.1 -0.05 



0.05 0.1 0.15 0.2 



(a) 



(b) 



Figure 6. (a) Contour lines of tlie analytic solution p{x, t) explicitly given in equation 



(26). (b) Probability distribution P{p) (24). Dashed lines display the scalings of 



equations (30) and (31). 



in equation (28), we obtain at leading order (dropping the bar over p) 



Pt{p) 



P' 



with a = infF(A) = 2. 

A 



(30) 



This scahng is clearly observed on figure 6(b), where equation ([t]) is integrated using 
the stationary environment a{x,t) = sinx. 

Note that a clear power-law scaling is also seen on figure |6jb at small values of 
the mass density. To explain this scaling, let us to consider a cell of length Ax located 
at Xq far away from a zero of a{x) and where the density reaches its minimum. At 
leading order the mass ejected from this cell between t and t + dt is proportional to 
J ~ p{x,t)Tdx<J^{x)\:i:=xo as dxP ~ at this point. Considering that d^cr'^ix) is constant 



near xq yields an exponential decay of the mass in the cell. Introducing in equation (24) 
a density of mass with an exponential decay p{x, t) ~ poc'^'o^ (with a^ an effective rate) 
directly leads to 

PAp) ~ -, (31) 

P 

which is in agreement with figure [6]b. 

As we will now see the situation is rather different in the non-stationary case where 

the motion of the zeros of the ejection rate limits the process of mass concentration. Now, 



the coefficients ar{t) and ai{t) of equation (13) are the Ornstein-Uhlenbeck processes 



defined in equations (^11| and (12). The PDF of p for different values of T is presented 



on figure [7| The time needed for the density to accumulate near the zero of the ejection 
rate a^ is of the order of F^^. Hence, for small values of F (i.e. Ku ^ 1), the fast 
temporal variations of a{x,t) do not leave enough time for the system to cumulate 
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Figure 7. Log-log plot of the PDF of the density for K = 1 and different values of 
r. Inset: Same figure in log-lin. 



mass. It is apparent in this case that most of mass is distributed around the mean 
mass (p) = 1 (see the F = 0.1 curve represented by orange pentagrams). However, 
for large but finite values of Ku, accumulation is fast enough and a resulting p~^^'^ 
scaling is observed. This scaling can be derived by considering that for large Ku the 
system rapidly relaxes to a quasi-stationary solution in a time of order l/F and stays 
there for a time of the order of the correlation time of a. Far enough from a zero, the 
environment behaves as cr{x) ~ Cx and its temporal evolution can be neglected. The 
density converges there to a quasi-stationary solution, which is such that 9^(o"^p) ~ 0, 
so that p ~ (T^^ oc 1/x^. Closer to the zero, the time variations of its location become 
faster than the concentration of mass and the density saturates to a finite value Pmax 
(see the left-hand side of figure [8|. The transition between these two regimes occurs at 




t 


Ax ^ 
^ s- -^ 


T 


Ijl 







Figure 8. Left: sketch of the density field in the vicinity of a zero of the ejection 
rate a^ located at x = 0; the density grows as x~'^ and saturates to Pmax at a distance 
Ax ^ Pmax of the zero. Right: sketch of the diffusive time evolution of a zero of a^, 
which has remained at a distance less than Ax for a time T before the reference time 
t = 0. 

— 1/2 

a distance Ax from the zero, which by continuity satisfies Ax ~ pmax . The length Ax 
is of the order of the distance travelled by the zero during a timescale equal to that of 
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mass concentration. Hence, if we assume that the zero diffuses, one has Aa; ~ F"^/^, 
so that the typical value of Pmax is of the order of F. The distribution of density has 
thus a crossover at p ~ F. When p ^ F, that is for values much below the plateau 
at Pmax, the behavior is dominated by the divergence of the quasi-stationary density 



profile. Introducing p(x) ~ 1/x^ in equation (24) straightforwardly leads to the algebraic 
behavior 

P{p) ~ p-3/2 for K p < F. (32) 

Note that this value of the exponent in the large-density intermediate tail is very robust. 
It has also been observed for larger numbers of modes and with a environment of the 
form o"(x, t) = vF cos (x — ct) (data not shown). 

The behavior at p ^ F of the PDF of density is related to the large fluctuations 
of Pmax and thus to the events when the zero does not move much during a time larger 
than F~^. More precisely, if at a fixed time we observe a large value of pmax, this means 

— 1/2 

that the zero has not moved by a distance larger than pmax during a time interval of 
length T larger than F~^ (see the right-hand side of figure p|. When the zeros diffuse, 
this probability relates to the first exit time distribution of a diffusive process. The 
probability density of the first time T at which the Wiener process exit the interval 
|x| < Aa; is ~ exp(— CT/Ax^) (see, e.g., [20]). We thus obtain 

/■oo 

Prob (pma. > p) ~ / e-^^^dT. (33) 



For p ^ F the leading-order behavior is given by T = F ^. This leads to the following 
exponential cutoff for the PDF of mass density 



P{p) ~e-^^" ''forp>F. (34) 

This exponential behavior is confirmed numerically as can been seen on the lin-log plot 
showed in the inset of figure [7| 

4-2. Non-smooth random environment 

We now turn to the case of a time-dependent non-smooth ejection rate a. Figure |9] 
presents the PDF of p for a number of runs with different values of F and h. Remark 
that for large F a power-law behavior is observed at small masses — see inset of figure 
9(b)[ These tails can be understood using a phenomenological argument similar to the 



one used in [T7|. Consider an extreme event leading to a very low density in a cell 
where only ejection of mass has taken place for some time. For such a configuration 
we expect an exponential decay of the density and thus a time of order T oc — F^^ Inp 
to reach a small density p (where ctq is an effective rate). The probability of having 
this configuration depends only on the properties of the environment. Let p < 1 be the 
probability of having such a configuration. As the environment de-correlates in a time 
order one, the number of independent realizations of the environment that is needed 
for the full process to take place during a time T, is oc T oc — F~^lnp. Therefore, the 
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(b) 

9. Log-log plot of the PDF of the density p for different values of h and (a) 
and (b) F = 10. 



probability of this complete event is p-'-^^^p/^^ where C is a proportionality constant 
that depends on h. We thus obtain that the probability of such an extreme event is 
proportional to p'^, with (3 = —C\np/T > 0. Note that this is actually a lower bound 
to the small mass tail, which is valid as far as p 7^ and C 7^ 0. Obtaining an analytic 
expression of the exponent /3 is a theoretically challenging problem that is beyond the 
scope of the present work. The exponent /3 can be obtained by fitting the left tail of 
the mass density PDF inside the range where the power law is observed. The fits are 



presented in dashed lines on the inset of figure 9(b) and the dependence of {3 on the 
Holder exponent h is displayed on figure 10(a) for T = 10. 




(a) 




(b) 



Figure 10. (a) exponent /? of the power-law fit at the small-mass tail of the PDF's 
on 9(b) i.e. for F = 10 and various h as labeled, (b) Variance {Sp'^) = {{p — (p))^) as 



function of h for different values of F. 
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Note that width of the PDF's on figure |9] clearly depends on the value of the Holder 
exponent h as it is apparent on figure [lOjb that represents the variance of the density 
{Sp"^) = {{p— (p))^)- Observe that both the exponent /3 and the variance 6p, do not have 
a monotonic behavior as a function of h. This can be qualitatively understood in the 
following way. When decreasing the value of h, the number of zeros of a{x, t) increases. 
However their spatial distribution also depends on the Holder exponent h. On the one 
hand, for h > 1/2 the space increments of a are positively correlated. This implies that 
there is no accumulation of zeros. For this configuration, the mass is transported from 
the non-vanishing zones of a{x,t) to the nearest zeros. Therefore decreasing h, more 
and more zeros are present, creating large masses and void zones and thus increasing 
fiuctuations. On the other hand, for h < 1/2 the covariance of increments is negative 
and there are some finite-size regions with a large number of zeros. cr(a;, t) vanishes many 
times in such zones, so that the diffusion is very weak and mass is trapped. This makes 
the transfer inefficient and reduces the probability of having extremely large or low mass 
concentrations. In simpler words, a large number of zeros increases fluctuations as long 
as they are not too dense. We expect then a change of behavior near oi h = 1/2, this is 



in agreement with flgure 10 a-b, especially for the largest value of T where the properties 
of the system are expected to depend more strongly on the ejection rate distribution of 
zeros. 

5. Scale invariance of the mass density field 

In this section we study the scaling properties of the density field. For that we introduce 
the coarse-grained mass density 

^ l'x+e/2 

Pe{x,t) = - p{y,t)dy. (35) 

^ Jx~e/2 

Note that by homogeneity we have (p) = (pe). The PDF of pe is computed using the 



definition (24) and it is displayed on figure 11 for different values of h and F = 1. The 



narrowest curves on each figure correspond to the largest (non trivial) scale L/2 and 
the curves monotonically increase their width when decreasing the scale down to the 
smallest scale of the runs. 

For spatially smooth environments, the coarse-grained density of mass becomes 
invariant when decreasing the scale, as it is apparent on figure 11 (a) [for h = 2. This 



indicates that the density field is spatially smooth and that the zeros are isolated. The 
collapse is faster for small densities. This is due to the accumulation of mass only in 
some very small clusters and the creation of large void zones. This large voids dominate 
the coarse-grained density statistics up to their typical size, that is for i < L/8, as 



seen from figure 11(a) On the contrary, the typical size 4iust of mass clusters is rather 
small. Averaging over scales larger than that this length scale reduces the largest mass 
fiuctuations. Indeed, the coarse-grained density cannot exceed 4iust/^ and this cutoff 



decreases with i, as can be observed on figure 11 (a 
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Figure 11. PDFs of the density p for F = 1, different h and various coarse-graining 
scales. The narrowest curve on each figure corresponds to the largest scale L/2 and 
the curves monotonically increase their width when decreasing the scale down to the 
smallest scale of the runs, i.e. L/2^^. 



(a) 





(b) 



Figure 12. . (a) Variance {Spj) — {{p^ — {p})'^) a-s a function of ^ for different values 
of h, as labeled. Inset: exponent ({h;T) of the variance (Spj) ~ ^-CC'-r) for F = 1. 
b) PDF of the coarse-grained density p£ for ft, = 1/2 and different values of F and £ 



showing a collapse at large masses when the scaling ( 36 ) is used. 



This argument does not applies for non-smooth environments where the distribution 
of zeros is highly non trivial. The scale invariance is then broken, as observed on figures 



11(b) and 11(c) This is also apparent on figure 12(a) which represents the variance 
^p'e = {{Pi ~ {p)Y) ^s a function of i for different values of h. For h > 1 and £ <^ L the 
variance 6pj does not depend on the scale £ (see the inset of the figure). However, one 
observes that for h < 1 and £ <^ L the variance presents a power-law scaling 5p1 ~ £~^ 
that suggests a self-similar behavior of the density. The exponent C,{h] V) is displayed 



on the inset of figure |ll(c^ 



The scaling invariance of the density field can be understood in the neighborhood 
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of the zeros of a. Let us assume that the ejection rate vanishes at a; = xq. Then, in 
the neighborhood of Xq, we have a{x,t)^ = \a{xo + x,t)p ~ x'^^. Supposing that the 
density behaves as a power-law in the vicinity of xq and rescahng space as x = Aa;, one 
can easily see that the rescaled density p = p(x) is a solution of 

^ = rA--^[a^(x,t)p(x,t)]. (36) 

Hence, we expect that if the coarse-grained density pe presents a self-similar property, 
then the PDF of pe obtained with a given value of F will coincide at large values with 
that of p\e corresponding to an ejection rate of amplitude FA^"^''. This scale invariance 



is apparent on figure 12 b where different PDF's of pe, with different values of F and 



such that Ti"^^ ^ = const are confronted for h = 1/2 and collapse at large densities. 

6. Conclusions 

To summarize, we have introduced a diffusion model based on a simple discrete mass 
ejection process where the ejection rates are random variables with temporal and 
spatial correlations. The model can be interpreted in terms of random walks in a 
time-dependent random environment. We have considered space-periodic environments 
where the temporal dependence of each Fourier mode is given by independent Ornstein- 
Uhlenbeck processes and both the amplitudes of the modes and their correlation times 
present some prescribed spatial scaling properties. This allowed us to consider smooth 
and non-smooth environments with fast and slow temporal dependence. No other 
assumptions were made on the environment and we expect such environments to display 
generic properties and to be representative of sufficiently general situations. 

The model was studied analytically and numerically. The corresponding particle 
dynamics is given by an Ito differential equation with a multiplicative noise and no 
drift. We observed that random trajectories diffuse at large times. Also, the probability 
distribution of their displacements tends to a Gaussian at large times and the deviations 
to this asymptotic behavior decrease as t~^. To take into account the fluctuations 
of mass due to the randomness of the environment we introduced and studied the 
probability distribution of the density of mass. We obtained some analytical results 
on the density of mass in the case of a stationary environment and showed that it 
displays a power-law tail at large masses with exponent —2. In the general case, we 
observed a competition between a trapping effect due to vanishing ejection rates of the 
random environment and the mixing due to its temporal dependence which leads to 
large fluctuations of the density of mass. These fluctuations were studied for both, 
smooth and non-smooth random environment. In the smooth case, we showed that the 
PDF has an intermediate power-law behavior, like in the stationary case but with an 
exponent —3/2, followed by an exponential cutoff. Finally, we studied the spatial scaling 
properties of the mass distribution by introducing a coarse-grained density fleld. For 
smooth random environments the coarse-grained density was found to be scale invariant. 
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We showed that at large masses, it possesses some scahng properties that depend on 
the coarse-graining size, the ejection rate typical amplitude, and the Holder exponent 
of the environment. 

The overall dynamics of the model that was introduced in this paper contains several 
space and time scales. Mass rapidly accumulates near the regions with a vanishing 
ejection rate and slowly moves following the diffusion of the zeros. Depending on the 
properties of the environment, there is also a clear separation of length scales: small mass 
fluctuations at small scales and large ones at the scale of the distance of two vanishing 
ejection rates. This scale separation strongly suggests trying to determine by standard 
homogenization techniques, a large-scale effective diffusion tensor and eventually an 
effective transport term. This matter is kept for future work. 

Extending our approach to dimensions higher than one is another possible direction 
for future. The ejection rate will then vanish on complicated sets that, depending on its 
regularity, might display fractal properties. Varying the Holder exponent of the random 
environment, its amplitude and the observation scales could then lead to a rather rich 
collection of different regimes. For instance, a particular attention would be needed to 
understand whether or not a power-law is still present in the mass density probability 
distribution. 

To conclude, let us stress again that the zeros of the environment play a crucial role 
on both the diffusive properties and the mass density statistics. They are responsible for 
mass accumulation but, at the same time, constitute barriers to transport. This duality 
implies that the fine statistical details of such systems crucially depend on the features 
of such zeros and in particular on the local structure of the ejection rate in their vicinity. 
In this work, we have decided to focus on ejection rates that can be written as the square 
of a generic Gaussian random function. However this choice cannot be always relevant 
and, for instance, it might be sometimes more appropriate to write the ejection rate 
rather as the exponential of a random function. This would change drastically most of 
the results reported in this paper. 
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